6 Amphibians
6.1 Calotriton asper
6.1.1 Retrieve data
species="Calotriton asper"
genus=strsplit(species, " ")[[1]][1] %>% tolower()
sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
rename(dataset=Dataset) %>%
filter(Species == species)
read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0121/DMB0121_counts.tsv.gz") %>%
rename(genome = 1)
genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0121/DMB0121_mag_info.tsv.gz") %>%
rename(genome = 1, length=mag_size)
genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0121/DMB0121_coverage.tsv.gz") %>%
rename(genome = 1)
download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0121/DMB0121.tree.gz", "data/DMB0121.tree.gz")
genome_tree <- read_tree("data/DMB0121.tree.gz")6.1.2 Filter data
#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
mutate(across(-1, ~ . * read_counts[[cur_column()]]))
# Transform to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))
# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))6.1.3 Community barplot
# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Stacked barplot
genome_counts %>%
mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE))) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
labs(y = "Relative abundance") +
facet_nested(. ~ Sample + Extraction, scales="free_x") +
guides(fill = guide_legend(ncol = 3)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
legend.position="none",
legend.title=element_blank())6.1.4 Alpha diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=0) %>%
t() %>%
as.data.frame() %>%
rename(richness=1) %>%
rownames_to_column(var="dataset")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1) %>%
t() %>%
as.data.frame() %>%
rename(neutral=1) %>%
rownames_to_column(var="dataset")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1,tree=genome_tree) %>%
t() %>%
as.data.frame() %>%
rename(phylogenetic=1) %>%
rownames_to_column(var="dataset")
# Merge alpha diversities
alpha_diversity <- richness %>%
full_join(neutral,by=join_by(dataset==dataset)) %>%
full_join(phylogenetic,by=join_by(dataset==dataset))
# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))
# Print alpha diversity
alpha_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
group_by(Extraction) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| Extraction | richness | neutral | phylogenetic |
|---|---|---|---|
| DREX | 107.0 | 51.75880 | 6.406283 |
| EHEX | 106.5 | 51.56757 | 6.311931 |
| ZYMO | 97.5 | 48.75158 | 6.219160 |
6.1.5 Beta diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C", out="pair")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C", out="pair")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")
# Merge beta diversities
beta_diversity <- richness %>%
full_join(neutral,by=c("first", "second")) %>%
full_join(phylogenetic,by=c("first", "second")) %>%
rename(richness=C.x, neutral=C.y, phylogenetic=C)
# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))
# Print beta diversities
beta_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
mutate(relationship = case_when(
Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
TRUE ~ NA_character_ # Handle other cases if needed
)) %>%
filter(relationship != is.na(relationship)) %>%
group_by(relationship) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| relationship | richness | neutral | phylogenetic |
|---|---|---|---|
| inter | 0.29497046 | 0.381599435 | 0.105560612 |
| intra | 0.03215577 | 0.007477988 | 0.001740429 |
6.1.6 Variance partitioning
richness <-genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
write_tsv(paste0("results/var_",genus,".tsv"))6.2 Lissotriton helveticus
6.2.1 Retrieve data
species="Lissotriton helveticus"
genus=strsplit(species, " ")[[1]][1] %>% tolower()
sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
rename(dataset=Dataset) %>%
filter(Species == species)
read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0124/DMB0124_counts.tsv.gz") %>%
rename(genome = 1)
genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0124/DMB0124_mag_info.tsv.gz") %>%
rename(genome = 1, length=mag_size)
genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0124/DMB0124_coverage.tsv.gz") %>%
rename(genome = 1)
download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0124/DMB0124.tree.gz", "data/DMB0124.tree.gz")
genome_tree <- read_tree("data/DMB0124.tree.gz")6.2.2 Filter data
#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
mutate(across(-1, ~ . * read_counts[[cur_column()]]))
# Transform to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))
# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))6.2.3 Community barplot
# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Stacked barplot
genome_counts %>%
mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE))) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
labs(y = "Relative abundance") +
facet_nested(. ~ Sample + Extraction, scales="free_x") +
guides(fill = guide_legend(ncol = 3)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
legend.position="none",
legend.title=element_blank())6.2.4 Alpha diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=0) %>%
t() %>%
as.data.frame() %>%
rename(richness=1) %>%
rownames_to_column(var="dataset")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1) %>%
t() %>%
as.data.frame() %>%
rename(neutral=1) %>%
rownames_to_column(var="dataset")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1,tree=genome_tree) %>%
t() %>%
as.data.frame() %>%
rename(phylogenetic=1) %>%
rownames_to_column(var="dataset")
# Merge alpha diversities
alpha_diversity <- richness %>%
full_join(neutral,by=join_by(dataset==dataset)) %>%
full_join(phylogenetic,by=join_by(dataset==dataset))
# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))
# Print alpha diversity
alpha_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
group_by(Extraction) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| Extraction | richness | neutral | phylogenetic |
|---|---|---|---|
| DREX | 22 | 12.05083 | 4.439034 |
| EHEX | 22 | 12.15999 | 4.408855 |
| ZYMO | 21 | 10.73360 | 3.336904 |
6.2.5 Beta diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C", out="pair")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C", out="pair")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")
# Merge beta diversities
beta_diversity <- richness %>%
full_join(neutral,by=c("first", "second")) %>%
full_join(phylogenetic,by=c("first", "second")) %>%
rename(richness=C.x, neutral=C.y, phylogenetic=C)
# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))
# Print beta diversities
beta_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
mutate(relationship = case_when(
Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
TRUE ~ NA_character_ # Handle other cases if needed
)) %>%
filter(relationship != is.na(relationship)) %>%
group_by(relationship) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| relationship | richness | neutral | phylogenetic |
|---|---|---|---|
| inter | 0.9090909 | 0.9487207 | 0.3681096 |
| intra | 0.2288066 | 0.1798833 | 0.1289869 |
6.2.6 Variance partitioning
richness <-genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
write_tsv(paste0("results/var_",genus,".tsv"))6.3 Salamandra atra
6.3.1 Retrieve data
species="Salamandra atra"
genus=strsplit(species, " ")[[1]][1] %>% tolower()
sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
rename(dataset=Dataset) %>%
filter(Species == species)
read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0129/DMB0129_counts.tsv.gz") %>%
rename(genome = 1)
genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0129/DMB0129_mag_info.tsv.gz") %>%
rename(genome = 1, length=mag_size)
genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0129/DMB0129_coverage.tsv.gz") %>%
rename(genome = 1)
download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0129/DMB0129.tree.gz", "data/DMB0129.tree.gz")
genome_tree <- read_tree("data/DMB0129.tree.gz")6.3.2 Filter data
#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
mutate(across(-1, ~ . * read_counts[[cur_column()]]))
# Transform to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))
# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))6.3.3 Community barplot
# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
select(colors) %>%
pull()
# Stacked barplot
genome_counts %>%
mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE))) %>% #apply TSS nornalisation
pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
scale_fill_manual(values=phylum_colors) +
labs(y = "Relative abundance") +
facet_nested(. ~ Sample + Extraction, scales="free_x") +
guides(fill = guide_legend(ncol = 3)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
axis.title.x = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
legend.position="none",
legend.title=element_blank())6.3.4 Alpha diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=0) %>%
t() %>%
as.data.frame() %>%
rename(richness=1) %>%
rownames_to_column(var="dataset")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1) %>%
t() %>%
as.data.frame() %>%
rename(neutral=1) %>%
rownames_to_column(var="dataset")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1,tree=genome_tree) %>%
t() %>%
as.data.frame() %>%
rename(phylogenetic=1) %>%
rownames_to_column(var="dataset")
# Merge alpha diversities
alpha_diversity <- richness %>%
full_join(neutral,by=join_by(dataset==dataset)) %>%
full_join(phylogenetic,by=join_by(dataset==dataset))
# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))
# Print alpha diversity
alpha_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
group_by(Extraction) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| Extraction | richness | neutral | phylogenetic |
|---|---|---|---|
| DREX | 70.0 | 35.33416 | 6.024485 |
| EHEX | 84.0 | 38.29817 | 6.183859 |
| ZYMO | 58.5 | 30.15278 | 5.653610 |
6.3.5 Beta diversity
#Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C", out="pair")
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C", out="pair")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")
# Merge beta diversities
beta_diversity <- richness %>%
full_join(neutral,by=c("first", "second")) %>%
full_join(phylogenetic,by=c("first", "second")) %>%
rename(richness=C.x, neutral=C.y, phylogenetic=C)
# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))
# Print beta diversities
beta_diversity %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
mutate(relationship = case_when(
Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
TRUE ~ NA_character_ # Handle other cases if needed
)) %>%
filter(relationship != is.na(relationship)) %>%
group_by(relationship) %>%
summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
tt()| relationship | richness | neutral | phylogenetic |
|---|---|---|---|
| inter | 0.3854701 | 0.3862921 | 0.125862853 |
| intra | 0.1422510 | 0.0315443 | 0.007370115 |
6.3.6 Variance partitioning
richness <-genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=0, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
neutral <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hillpair(.,q=1, tree=genome_tree, metric="C") %>%
as.dist() %>%
adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
broom::tidy() %>%
filter(term == "Extraction") %>%
select(R2) %>% pull()
tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
write_tsv(paste0("results/var_",genus,".tsv"))